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ABSTRACT 

We study the combined effects of relaxation, tidal heating and binary heating on 
globular cluster evolution, exploring the physical consequences of external effects and 
examining evolutionary trends in the Milky Way population. Our analysis demonstrates 
that heating on circular and low-eccentricity orbits can dominate cluster evolution. The 
results also predict rapid evolution on eccentric orbits either due to strong relaxation 
caused by the high densities needed for tidal limitation or due to efficient bulge shocking 
of low density clusters. 

The combination of effects leads to strong evolution of the population as a whole. 
For example, within the solar circle, tidally-limited W^Mq clusters lose at least 40% of 
their mass in 10 Gyr. At high eccentricity most of these clusters evaporate completely. 
Bulge shocking disrupts clusters within 40 kpc which have less than 80% of their mass 
within their pericentric inner Lagrange point. Our results are consistent with sugges- 
tions that the shape of the cluster luminosity function results from evaporation and 
disruption of low mass clusters; they further predict that the net velocity dispersion of 
the cluster system in the inner Galaxy has decreased with time. Preliminary constraints 
on formation models are also discussed. We conclude that the observed cluster system 
has largely been shaped by dynamical selection. 
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^ 1 INTRODUCTION 



Many studies of globular cluster evolution have focused on 
internal mechanisms which drive evolution. This work has 
produced a clear picture in which initial stellar evolution 
causes significant mass loss from a nascent cluster (e.g. Cher- 
nofF & Weinberg 1990); two-body relaxation leads to mass 
segregation (e.g. Inagaki & Saslaw 1985) and core collapse 
in surviving clusters (e.g. Cohn 1980); binary heating halts 
collapse (e.g. Lee et al. 1991); and the cluster continuously 
loses mass due to the escape of stars, eventually undergoing 
complete evaporation (e.g. Lee & Goodman 1995). 

It is also recognized that the Galaxy influences clus- 
ter evolution. The time-dependent tidal field heats clus- 
ters and tidal limitation aids in the removal of escaping 
stars. Previous investigations have considered disk shocking, 
bulge shocking and tidal limitation, concluding that each 
will play a role, particularly in the inner Galaxy (e.g. Spitzer 
& Chevaher 1973; Chernoff & Shapiro 1987; Aguilar, Hut & 
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Ostriker 1988; Weinberg 1994). In addition, recent observa- 
tional studies showing correlations of cluster properties with 
Galactocentric position indicate the measurable infiuence of 
the Galaxy (e.g Chernofi' & Djorgovski 1989; Djorgovski et 
al 1993; Djorgovski & Meylan 1994). 

The principal tool used in studies of cluster evolution 
over the last decade-and-a-half has been direct solution of 
the Fokker-Planck equation (Cohn 1979). However, most 
of these calculations have excluded external effects. Re- 
cently, using time-dependent perturbation theory to investi- 
gate disk shocking, Weinberg (1994) demonstrated that res- 
onant forcing can perturb internal stellar trajectories beyond 
the limit expected from adiabatic invariance. This indicates 
that the Galaxy plays a greater role in cluster evolution than 
previously thought and motivates new studies of cluster evo- 
lution which combine internal and external effects. 

The importance of external heating requires us to re- 
examine the current picture of internally-driven evolution. 
In particular, external effects will infiuence the collapse 
rates, evaporation times and general physical properties de- 
rived in previous calculations. The present work compares 
this behavior with and without heating over a wide range of 
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cluster properties to present a revised view. This study also 
examines the survival and disruption characteristics of clus- 
ters on a range of Galactic orbits to shed light on the initial 
conditions of the cluster system. The results demonstrate 
that evolution does indeed depend strongly on position and 
orbit, further implying that observed cluster properties have 
been largely determined through dynamics. 

Our study rests on a linear theory of external heating- 
based on Weinberg's (1994) treatment of disk shocking- 
which we include in numerical solutions of the Fokker-Planck 
equation. Nearly all previous work has emphasized impulsive 
shock heating due to a single passage through the disk or 
bulge. The work presented here describes resonant heating 
due to the time-varying tidal field encountered on periodic 
orbits of the cluster in the Galaxy- an effect we refer to as 
orbit heating. In this context, shock heating is seen to result 
from the broad resonances caused by an impulsively applied 
external force. 

Although our treatment of external heating can include 
the influence of any component of the Galactic potential, 
here we consider only the spheroid in order to allow precise 
definition of the physical behavior and preliminary descrip- 
tion of the evolutionary trends. The present study includes 
heating on cluster orbits in the isothermal sphere and is used 
to study cluster evolution from initial King model states to 
the point of complete evaporation on a range of orbits in the 
Galaxy. Our conclusions, therefore, place only lower limits 
on the overall rate of cluster evolution but are significant 
nonetheless. 

The plan of the paper is as follows. We derive the lin- 
ear theory of external heating in §^ and discuss its physical 
interpretation in In §^ the numerical implementation is 
described. In §^ we present the results of our study of clus- 
ter evolution under the combined influence of internal and 
external effects. Finally, in §^ we discuss the implications of 
the results for the Milky Way globulars. Readers concerned 
primarily with the effects of heating and its evolutionary 
consequences may skip M without loss of continuity. 



2 DERIVATION OF EXTERNAL HEATING 
RATE 

The physics behind the perturbation theory discussed be- 
low can be summarized as follows. Each orbit in the cluster 
acts like a pendulum with two-degrees of freedom (cf. Bin- 
ney & Tremaine 1987, Chap. 3). The time-dependent tidal 
fleld can drive the pendula at a discrete or continuous spec- 
trum of frequencies depending on whether the perturbation 
is quasi-periodic or aperiodic, respectively. Because the tem- 
poral variation discussed here is caused by the cluster's orbit 
in the spherical component of the Galaxy, the spectrum is 
discrete. For disk shocking described by Weinberg (1994), 
the spectrum is continuous. In both cases, the energy of 
each orbit changes as it passes through the resonance. The 
accumulated effect of all possible resonances on all orbits, 
drives the secular evolution of the equilbrium distribution 
function (DF). The expressions given below are valid for 
both periodic and aperiodic cases. 

We compute the evolution by expanding the Boltzmann 
equation to flrst order and solving for the perturbed dis- 
tribution function (neglecting self-gravity). The flrst-order 



change phase mixes but second order energy input leads to 
an induced phase space flux which helps drive cluster evo- 
lution. N-body comparisons shown in Appendix ^ indicate 
that the self-gravity of the tidally-induced wake has negligi- 
ble effect for cases of interest here. 

We use a locally inertial reference frame which is cen- 
tered on the cluster and has axes fixed in space (see Ap- 
pendix ^ for derivation). The unperturbed Hamiltonian 
is therefore completely separable, implying the existence 
of action-angle variables. This frame allows the internal 
dynamics to be defined in accordance with the standard 
Fokker-Planck technique (e.g. Cohn 1979) which uses an 
energy-space DF f(E) and depends on the adiabatic invari- 
ance of the radial action. Within this framework, we derive 
a version of the formalism presented by Weinberg (1994) 
which describes heating of globular clusters on arbitrary or- 
bits in external potentials. 



2.1 Perturbed distribution function 

The linearized Boltzmann equation is a linear partial differ- 
ential equation in seven variables. Using action-angle vari- 
ables, we can separate the equation and employ standard 
DFs constructed according to Jeans' theorem (Binney & 
Tremaine 1987). The explicit form of the linearized Boltz- 
mann equation is 



dfi ^ dfi dHo dfo dH, _ ^ 

dt 9w dl dl 9w ' 



(1) 



where w is the vector of angles, and I are the conjugate ac- 
tions. The quantities /o and Hq depend on the actions alone. 
The small variation in Galactic potential over a typical clus- 
ter size allows quadratic expansion of the tidal field (see 
Appendix^ for details). We may thus define Hi — u{r)g{t) 
and expand in a Fourier series in action-angle variables (e.g. 
Tremaine & Weinberg 1984). Each term /n in the Fourier 
series is the solution of the following differential equation: 



^-Kd.f^)fH=il.§U,(I)g(t)^il.§Hn, 

where fl = dH/dl and 
1 



(27r)3 



/ \ ~il-w i3 

u(rje d w. 



(2) 



(3) 



The inhomogeneous equation may be solved using a Green's 
function (e.g. Birkhoff & Rota 1962, p. 39) to give the time- 
dependence for each coefficient of the perturbed DF 



/ii = il-^U,(I)e- 



dt e g(t ) 



(4) 



where we have assumed that the perturbation begins at time 
to. 



2.2 Heating rate 

The rate of change in energy arising from the perturbation 
follows from Hamilton's equations (Weinberg 1994). The to- 
tal phase-averaged change in energy can be written as 



(E) 



dtJ2 (il-f^)Hi-ifii. 



(5) 
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Substituting for /n from equation yields 

(£) = -4.=' 5](i.f2)(i.^)|u,r / dtV'--'g(t'; 



■ (6) 



which represents the heat input due to the perturbation dur- 
ing an interval = t — to- This expression is valid for finite- 
duration, aperiodic perturbations such as disk passage as 
well as periodic perturbations which arise on regular orbits 
in the Galaxy. In particular, Weinberg's (1994) results for 
disk shocking are obtained from equation ^ by substitut- 
ing the tidal amplitude appropriate to the disk profile for 
g{t') and integrating over the interval (—00,00) assuming a 
linear trajectory. 

For periodic perturbations it is more suitable to derive 
the asymptotic heating rate (e.g. Landau & Lifschitz 1965, 
p. 151). We first expand the tidal amplitude in a Fourier 



?(t) ^ ^ an 



(7) 



n— ~ OC 



and substitute into equation (^. Taking the limit t — > 00 
and assuming the onset of the perturbation at to = 0, we 
obtain 



n — — OC 



Integrating (E) over inclination and angular momentum, 
we obtain the change in energy which defines the induced 
change in the distribution, given by a one-dimensional con- 
tinuity equation in energy space (appropriate to the ID 
Fokker-Planck formulation employed below; see Appendix 
^ for derivation) : 



21 
dt 



d 



16tt^P{E) dE 



{{(E))}, 



(9) 



where P{E) is the phase space volume. This is called the 
equation of quasilinear diffusion in the plasma literature 
(e.g. Stix 1992). The term quasilinear refers to the propor- 
tionality of the heating rate to the squared amplitudes of 
the linear modes. The linear modes arise from the resonant 
forcing of stellar orbits by a periodic perturbation. The com- 
petition between two-body relaxation and this externally in- 
duced phase space fiux can strongly infiuence globular clus- 
ter evolution, as we will show below. 



(11) in Weinberg (1994), we can write the perturbation as a 
series in action-angle variables; 



00 

Hi = ij^oW J2 {^^/4^^ooo(/3)<5^,o<5^30 



^ ^3 

1— — 00 



-e^"^/^v2,„-2(/3)^.3,-2}x; 



e 



where 



X',' = — 



dwie 



(11) 



(12) 



/3 is the inclination of the orbital plane and Vu^i^iP) is a 
rotation matrix (e.g. Tremaine & Weinberg 1984) . The angle 
'ip—W2 is the difference between the mean azimuthal angle W2 
and the azimuthal angle in the orbital plane. We substitute 
this expansion into equation (^) and average over inclination 
and angular momentum to derive the heating rate 



OC „ 



a„ 1^5(1 ■ f2 — nuj) 



n— — 00 

OC 

n — — OC 



where 



an = - 



b„ = 



P/2 

■P/2 
P/2 

■P/2 



dtQ,l{t)e 



^S{l-n-nuj) }, (13) 

(14) 
(15) 



and P is the period of the cluster orbit. 

For an isotropic DF, 1 • dio/dl = (1 ■ n)dfo/dE. We also 
explore the effect of anisotropy using Merritt-Osipkov mod- 
els (e.g. Binney & Tremaine 1987). The distribution func- 
tion takes the form /o(I) = f(Q), where Q ^ E± /2rl, 
1-afo/ai = dfo/dQ{\ ■ n =f li^^iJMr^ ± hl/rl). The 
anisotropy increases with decreasing anisotropy radius, ra- 



2.3 Heating rate in isothermal sphere 

Below we will need the heating rate for a cluster orbiting 
in the isothermal sphere. For most galaxies, the small varia- 
tion in potential over a typical cluster size allows quadratic 
expansion of the tidal field. Therefore, the perturbing Hamil- 
tonian is: 

Hi = if7o(i) [-cos2e(t)a;^ 2 sin 20 {t)xy 

- cos 20 {t)y'^ + z'^], (10) 

where flo{t) = Vo/R{t) is the angular rotation speed at the 
orbital radius at time t and Q{t) = J^^dt'Q is the instan- 
taneous azimuthal angle of orbit. Using equations (10) and 



3 DISCUSSION OF PHYSICAL MECHANISM 

A cluster orbiting in the Galaxy feels a time-dependent tidal 
field. A typical orbit is periodic and introduces a periodic 
external force on orbits of cluster stars. As described in §^ 
resonant heating occurs when the periods of stellar orbit 
and external force coincide, leading to repeated acceleration 
and increase in the energy of individual orbits. Integrated 
over many periods, the energy gained by the orbit increases 
linearly with time (c.f. eq. ^) . Energy absorption eventually 
leads to the evolution of individual orbits, (see Appendix O 
for discussion and numerical implementation of finite dura- 
tion resonances). This in turn drives the secular evolution 
of the cluster potential. 
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Orbits can either gain or lose energy to the tidal field 
depending on the particular resonance. For example, in disk 
galaxies with flat rotation curves it is well-known that the 
inner Lindblad resonance loses energy to a perturbation 
while an outer Lindblad resonance gains energy. However, 
for isotropic distribution functions with df /dE > 0, the per- 
turbation always heats the system on average though some 
localized regions of phase space may lose energy. 

Non-resonant interaction has no net effect on an orbit. 
Successive maxima in the external force tend to accelerate 
and decelerate the star equally, leading to asymptotic can- 
cellation as long as the initial transients remain linear (i.e. 
do not alter the intrinsic frequency of the star with an ini- 
tial jolt). Over short times, non- resonant heating does occur 
because the time duration is insufficient for complete can- 
cellation to occur. 

Non-linear transient or impulsive heating leads to rapid 
change in orbital energies as a rapidly applied force 'kicks' 
a star regardless of its orbital frequency. However, the stan- 
dard impulse approximation, when used to describe a peri- 
odic perturbation, ignores the long-term decay of transient 
energy in the linear limit as well as the linear growth in en- 
ergy at the resonances. For most cases of interest, heating 
rates are in the linear limit, implying that external influence 
depends primarily on secular transfer of energy through or- 
bital resonances. 

To illustrate the behavior of transients and transient 
decay. Figure |l| compares the exact time-dependent energy 
input given by equation with the energy input defined by 
the asymptotic heating rate equation (^). Transients decay 
rapidly at low energy and more slowly at high energy. Empir- 
ically, we find that two to three Galactic orbital periods are 
required before the asymptotic limit is effectively reached. 
This treatment therefore adequately describes all but the 
outermost halo clusters for which initial transients may still 
be important. The comparisons of perturbation theory with 
N-body simulation shown in Appendix ^ demonstrate the 
validity of the approach. 

The magnitude of the heating rate is determined by the 
cluster profile, density and orbit. The profile and density de- 
fine the distribution of internal orbital frequencies and the 
cluster orbit defines the external forcing frequencies and am- 
plitude. For a cluster of fixed profile and mass, the density is 
determined by the tidal radius. Individual clusters may not 
be tidally limited due to initial conditions or heating-driven 
expansion. Therefore we use the function M{xp) to param- 
eterize the fraction of the total cluster mass enclosed within 
the instantaneous pericentric inner Lagrange point, Xp. This 
function depends on the profile and the ratio of cluster mean 
density to the mean density required by tidal limitation. 



A tidally-limited cluster has limiting radius, Rc 



and 



therefore M{xp) — 1, while a tidally-unlimited has Rc > Xp 
and therefore M{xp) < 1. Heating rates for a given orbit 
increase as M{xp) decreases. 

The perturbing potential in the logarithmic sphere, 
equation ([lo|), heats clusters on orbits of any eccentric- 
ity. The tide transfers energy and angular momentum to 
the cluster through the resonances, which unbinds stars. 
On circular orbits, the tidal field creates a triaxial pertur- 
bation of constant amplitude proportional to f^Q rotating 
with fixed pattern speed f2o- On eccentric orbits, conserva- 
tion of center-of-mass angular momentum introduces time- 
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Figure 1. The mean change in energy as a function of internal 
orbital energy in a Wq = 5 King model due to heating on an 
eccentric k = 0.3 (e 0.7) orbit after one orbital period. Com- 
parison of simulation (histogram), exact time-dependent pertur- 
bation theory (equation 6, solid) and heat input calculated from 
asymptotic heating rate (equation 8, dotted) shows that initial 
transients decay strongly at low energy while impulsive energy 
change persists at high energy. Horizontal dotted line indicates 
the level of accuracy in the simulation. 



dependent amplitude and rotation rate. This produces more 
resonances. Tidal torquing can also induce a net spin. 

The rate of external heating is also infiuenced by our 
choice of equilibrium phase space distributions. For exam- 
ple, according to Jeans' theorem, one can define equilibria 
in the rotating frame of a circular cluster orbit using the 
Jacobi constant, Ej (e.g. Heggie & Ramamani 1994). By 
transforming to the frame in which the perturbation is time- 
independent, we remove the resonances from the problem. 
We can therefore choose a bound distribution of orbits in 
Ej using the limiting zero-velocity surface, so no stars are 
lost and the cluster experiences no net tidal heating, al- 
though inertial energies and angular momenta are not con- 
served. Using f{E) instead of f{Ej) leads to heating in the 
analogous case because we cannot choose orbits which are 
strictly bound. In any case, a real cluster cannot reach true 
equilibrium because it is bound and therefore undergoes re- 
laxation. In fact, as is shown below, it is typically a compe- 
tition between external heating and relaxation due to strong 
resonances with diffused core stars that strongly influences 
cluster evolution. 



4 FOKKER-PLANCK CALCULATIONS 

To determine the influence of external heating on cluster 
evolution, we conduct a series of Fokker-Planck calcula- 
tions which begin with King model initial conditions and 
run through core collapse to complete evaporation. Relax- 
ation is computed using the multi-mass code of Chernoff & 
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Weinberg (1990) which solves Henon's (1961) orbit-averaged 
Fokker-Planck equation. Core heating is included in the form 
described by Lee et al (1991) with a time step that supresses 
stochastic core osillations. Irnplementation of external heat- 
ing is detailed in Appendix pi The comparisons shown in 
Appendix |e| are used to test the implementation. 

Each physical process depends on the input model pa- 
rameters listed in Table |^. The total mass is denoted by 
Mc and the concentration by Wo- Orbits in the isothermal 
sphere are defined by their energy E and angular momentum 
J. In place of absolute angular momentum J, we use the rel- 
ative angular momentum k = J / Jmax{E), where Jmax{E) 
is the angular momentum of a circular orbit with energy 
E. The value Hi = Q denotes a radial orbit and the value 
K = 1 denotes a circular orbit. The apocentric, pericentric 
and mean orbital radii are denoted Ra,Rp, Rav, respectively. 
We represent the Galactic potential as a singular isothermal 
sphere with rotation velocity vq — 220 km/s. 

We consider a range of initial values for M{xp). If the 
young, rich LMC clusters are representative of young globu- 
lar clusters, M{xp) may be significantly smaller than unity 
initially (E lson, Fall & Freeman 1987). Furthermore, as dis- 
cussed in §5.4.2, formation scenarios can imply varying de- 



grees of tidal truncation for an individual cluster depending 
on the local conditions under which it forms and the orbit 
on which it travels. 

The distribution of stellar masses in the cluster is given 
by a power-law mass spectrum, dN/dM oc m~°, with upper 
and lower mass limits mi and ruu, respectively. Fiduci al v al 
2.35 



5.1 



to 



ues a = 2. at), mi =0.1 and m„ = 2.0 are adopted in 
represent the cluster mass spectrum following the period of 
strong stellar evolution when relaxation, tidal heating and 
binary heating dominate cluster evolution. The importance 
of stellar evolution diminishes after ~ 1 Gyr for a = 2.35 and 
mi =0.1 which corresponds to the main sequence lifetime 
for a 2Mq A-st ar. The effect of changing the mass spectrum 
is explored in §5.2. 



5 RESULTS 

5.1 Orbital heating and bulge shocking 

Because heating rates depend on cluster profile, tidal trun- 
cation and orbit, comparisons in different physical regimes 
are needed to demonstrate the primary influences of heating 
on cluster evolution. We choose four specific examples listed 
in Table H 

Example 1 compares the relative strengths of heating 
rates on different orbits using physically identical clusters, 
each of which is tidally-limited at its orbital pericenter. In 
this case, because the average tidal strength is largest on 
circular orbits, heating is also strongest on circular orbits 
and decreases with eccentricity (Figure ^). 

To investigate the effect of heating on long-term evo- 
lution, we compare evaporation times, tev, for tidally lim- 
ited clusters of different mass, concentration and k. Table ^ 
shows tev normalized by the circular, IQ^Mq Wo ^J) case 



(arbitrary scaling to physical units is provided in §5.4.1). In 
these comparisons, clusters of a particular mass and con- 
centration are identical and clusters of differing mass and 
concentration possess the mean density required for tidal 
limitation. 
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Figure 2. Example 1: comparison of heating rates in identical 
tidally-limited Wq = 5 clusters on different orbits. Values of K 
are indicated to the right of each curve. Heating on circular or- 
bits dominates because the average tidal strength is highest, de- 
creasing with eccentricity because average tidal amplitude drops 
monotonically. Heating rates in circular and k = 0.3 case differ 
by about 2 orders of magnitude. In the circular case, orbits near 
the tidal boundary gain ~ 10% in energy in an orbital time. 



For identical clusters, tev decreases monotonically with 
K, reaching a minimum for circular orbits. Evaporation times 
can decrease by a factor of two in circular cases when tidal 
heating is included. The relative evaporation times refiect 
the relative strength of heating rates shown in Figure ^ 
Heating accelerates evolution because external forcing effi- 
ciently torques and expels high-energy core stars on radial 
orbits, as noted by Oh & Lin (1992) in N-body calculations. 
This reduces the local relaxation time in the core, enhances 
relaxation rates and causes rapid evaporation. Spitzer & 
Chevalier (1973) noted this effect in certain regimes of disk 
shocking, interpreting it as an increase in the core-halo tem- 
perature gradient (see also Chernoff & Shapiro 1987, Wein- 
berg 1994). For the highest eccentricities, tev is only slightly 
shorter than with no heating, demonstrating the insignifi- 
cance of high-eccentricity heating in tidally-limited clusters. 

In many cases, evaporation time does not vary strongly 
with concentration for the same orbit and mass, indicating 
that overall mass loss rates are insensitive to initial concen- 
tration. In the exceptional k = 0.9 and 1.0, Wo = 3, K^Mq 
cases, heating causes rapid disruption because these clusters 
have low binding energy and long relaxation times and are 
easily torn apart by the tide. 

While Example 1 compares heating rates as a function 
of eccentricity in identical clusters, the orbits occupy dif- 
ferent regions of the Galaxy (c.f. Table ^). In Example 2, 
we consider clusters in similar regions by comparing tidally- 
limited clusters on orbits of equal mean radius. Because they 
are tidally truncated, these clusters still undergo the same 
rate of heating relative to internal energies shown in Figure 
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Table 1. Processes and parameter dependences 



Process 



Parameters 

a mi rriu Wq M{xp) E k 



Relaxation \J \J \J \/ \/ \/ 

External Heating ^ \/ \/ \/ 

Core Heating V V \/ \/ 



Table 2. Example scenarios for a 10^ cluster 



Example 


K 


Ra (kpc) 


Rp (kpc) 


P (lOOMyr) 


rt (pc) 


Xp (pc) 


M{xp) 




1 


1.0 


8.5 


8.5 


2.5 


70 


70 


1.0 


5.0 




0.9 


11.2 


9.2 


2.1 


70 


70 


1.0 


5.0 




0.7 


19.6 


10.3 


3.1 


70 


70 


1.0 


5.0 




0.5 


37.8 


11.8 


5.3 


70 


70 


1.0 


5.0 




0.3 


89.4 


13.7 


11.3 


70 


70 


1.0 


5.0 


2 


1.0 


8.5 


8.5 


2.5 


70 


70 


1.0 


5.0 




0.9 


9.4 


7.6 


1.8 


63.5 


63.5 


1.0 


4.3 




0.7 


10.9 


5.7 


1.7 


48.5 


48.5 


1.0 


2.9 




0.5 


11.9 


3.7 


1.7 


33.2 


33.2 


1.0 


1.6 




0.3 


12.4 


1.9 


1.6 


19.3 


19.3 


1.0 


0.7 


3 


1.0 


8.5 


8.5 


2.5 


70 


70 


1.0 


5.0 




0.9 


9.4 


7.6 


1.8 


70 


61.4 


1.0 


5.0 




0.7 


10.9 


5.7 


1.7 


70 


47.6 


0.99 


5.0 




0.5 


11.9 


3.7 


1.7 


70 


31.5 


0.92 


5.0 




0.3 


12.4 


1.9 


1.6 


70 


16.4 


0.63 


5.0 


4 


0.3 


15.0 


2.3 


1.9 


22.9 


22.9 


1.0 


0.9 




0.3 


15.0 


2.3 


1.9 


41.3 


23.5 


0.95 


2.3 




0.3 


15.0 


2.3 


1.9 


48.8 


21.1 


0.9 


2.9 




0.3 


15.0 


2.3 


1.9 


61.1 


20.3 


0.8 


4.1 




0.3 


15.0 


2.3 


1.9 


73.2 


19.4 


0.7 


5.4 



^ However, cluster densities vary due to differences in or- 
bital angular frequencies. For fixed cluster mass, this implies 
that tidal radii will vary. 

The tidal radius rt decreases with the increased peri- 
galactic angular frequency at higher eccentricity. This in- 
creases the mean density and decreases the dynamical time 
tdyn, producing shorter relaxation times, larger evaporation 
rates and, as a result, shorter lifetimes as compared to Ex- 
ample 1. A cluster with k = 0.3 in Example 2 will evaporate 
in 1/7 the time of a cluster with k = 0.3 in Example 1 and 
1/5 the time of a cluster with k = 1.0 (the circular case). 

Since cluster orbits are generally unknown, the degree 
of tidal truncation at pericenter cannot be directly inferred. 
So, in Example 3, we assume that an observed cluster lies 
at its average orbital radius for a range of eccentricity and 
is tidally limited for zero eccentricity. The mass within the 
pericentric inner Lagrange point M{xp) can be substantially 
less than unity on eccentric orbits (Table ^). This leads to 
stronger heating than found on the same orbits in Example 
1 (see Figure 1^. For k = 0.7 the heating rate is much larger 
than in Example 1 even though only small amounts of mass 
overlie Xp. For k = 0.3 strong impulsive heating or bulge 
shocking (e.g. Aguilar et al 1988) occurs due to the increase 
in tidal amplitude. 



Example 4 shows the dependence of heating rates on 
degree of tidal truncation for a fixed k = 0.3 orbit. Figure ^ 
illustrates the dependence of heating on both k and M(xp): 
significant heating will occur on orbital timescales for k = 
0.3 and M{xp) < 0.9. Strong heating for k > 0.3 will also 
occur because these orbits have larger heating rates for the 
same value of M{xp). Table | shows the variation in cluster 
size and dynamical time with tidal truncation, indicating 
the corresponding variation in mean density. 

The evolutionary consequences of the heating rates in 
Example 4 are shown in Table ^ Clusters of different mass 
and concentration have equal M{xp) on the same orbit. 
Weakly bound clusters disrupt more easily because the res- 
onances occur more deeply within the system. Survival of 
Wo = 3 clusters decreases strongly with M{xp). Conversely, 
for small reduction in M{xp) , survival of Wo = 5 and Wo = 7 
clusters is enhanced as increased heating is offset by dimin- 
ished relaxation. For Wo = 5, maximum enhancement oc- 
curs at M(xp) ~ 0.95. For Wo = 7, higher binding energies 
lead to even longer lifetimes for more severe truncations. 
Further reductions in M{xp) eventually lead to rapid dis- 
ruption due to strong tidal shocking. 

These results define a rough criterion for bulge shock- 
ing: for Wo < 5 and k > 0.3, bulge shocking will occur for 
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Table 3. Evaporation times t, 



Wo = 3 






K 






Me (Mo) 


1.0 


0.9 


0.6 


0.3 


nh" 


1.0 X 10^ 
1.0 X 10^ 


0.65 
0.94 


0.95 
3.69 


1.27 
9.61 


1.29 
10.9 


I. 30 

II. 2 



Wo = 5 






K 






Me (Mo) 


1.0 


0.9 


0.6 


0.3 


nh 


1.0 X 10^ 


1.00 


1.15 


1.34 


1.37 


1.38 


5.0 X 10^ 


3.20 


4.25 


5.64 


6.10 


6.20 


1.0 X 10« 


5.20 


7.01 


10.31 


11.43 


11.82 



Wo = 7 






K 






Me (Mo) 


1.0 


0.9 


0.6 


0.3 


nh 


1.0 X 10^ 


1.12 


1.27 


1.40 


1.41 


1.42 


5.0 X 10^ 


3.97 


4.70 


5.95 


6.25 


6.37 


1.0 X 10^ 


6.35 


8.29 


11.06 


11.96 


12.13 



" nh denotes no heating 




Figure 3. Heating rates for Example 3. Tidally-limited circular 
case (dashed line) is plotted for reference. For k = 0.7, heating 
increases strongly in mildly tidally-unlimited cluster compared to 
Example 1. For k = 0.3, distribution with E > —0.3 undergoes 
strong shocking on orbital time scale. 



M{xp) < 0.9. Disruption for fixed M{xp) and k also impUes 
disruption for larger k because heating rates increase with 
K. For K < 0.3, bulge shocking requires even smaller M{xp) 
to cause disruption. 

This series of comparisons establishes three impor- 
tant aspects of tidal effects on different orbits in a spher- 
ical potential: 1) low-eccentricity and circular orbit heating 
for tidally-limited clusters strongly accelerate evolution; 2) 



Figure 4. Heating rates for Example 4. Shocking develops slowly 
as M{xp) decreases. For M{xp) ~ 0.95 heating of the tail is 
slightly stronger than in the tidally-limited circular case (dashed 
line). Heating at low energies is substantially less. Strong impul- 
sive heating or bulge shocking of the tail of the distribution will 
occur for M{xp) < 0.9. 



high-eccentricity heating has little effect in tidally-limited 
cases but the high mean density found for typical orbital 
radii in the Galaxy leads to short relaxation and evapora- 
tion times; 3) high-eccentricity heating, or bulge shocking, 
becomes important when clusters are tidally-unlimited, al- 
though the exact effect depends on M{xp), k, Wo and Mc- 
Finally, an important consequence of strong tidal heat- 
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Table 4. Bulge shocking evaporation times (Example 4) 



Wo = 3 






M{xp) 








1.0 


0.95 


0.9 


0.8 


0.7 


1.0 X 10^5 


1.3 


1.3 


1.1 


0.4 


0.2 


1.0 X 10^ 


10.9 


2.4 


0.8 


0.4 


0.4 


Wo = 5 






M{xp) 






M ( Mr-, ) 


1.0 


0.95 


0.9 


0.8 


0.7 


1.0 X 10^ 


1.4 


2.4 


2.1 


1.0 


0.5 


1.0 X 10^ 


11.4 


12.1 


5.2 


1.5 


0.5 


Wo = 7 






M{xp) 






Mc (Mq) 


1.0 


0.95 


0.9 


0.8 


0.7 


1.0 X 10^ 


1.4 


2.7 


2.4 


2.3 


2.0 


1.0 X lO'' 


12.0 


16.0 


23.4 


13.8 


7.5 



ing is suppression of the gravothermal instability. Although 
this may cause expansion and disruption, relaxation slows 
the expansion and can still produce mass segregation (Fig- 
ure^). Therefore, mass segregation does not necessarily im- 
ply core collapse, a possibility that does not arise when ne- 
glecting external heating (e.g. Chernoff & Weinberg 1990; 
Drukier, Richer & Fahlman 1992). Observed clusters with 
King profiles and strong mass segregation (such as M71) 
may reflect the influence of strong tidal effects. 

5.2 Influence of Mass Spectrum 

The mass spectrum controls the rate of relaxation and inter- 
play with external heating. Clusters with steep mass spectra 
or a narrow range of low mass stars have lower relaxation 
rates than do clusters with shallow mass spectra or a wide 
range in stellar mass (e.g. Chernoff & Weinberg 1990). Here 
we examine the competition between external heating and 
relaxation over a range in a and in unhealed and circu- 
larly heated tidally-limited clusters. 

Circular heating reduces evaporation times over a range 
m a (Table |). Roughly a factor of three reduction can oc- 
cur for a = 3.35. Differences between heated and unhealed 
clusters increase with a because the slower relaxation rates 
at high a are more readily enhanced. In addition, for fixed 
mass and concentration, heating reduces differences in evap- 
oration time which depend on a. 

Heating also reduces core collapse times tec up to 33% 
(Table |) and masses remaining at core collapse up to a 
factor of two (Table P). High concentration clusters maintain 
the same core collapse times in all cases but show decreased 
mass at core collapse. 

The non-monotonic behavior of core collapse time with 
spectral index was also found by Inagaki (1985 Table II) in 
Plummer law initial profile and Chernoff & Weinberg (1990 
Table 4) in King models. This indicates a complex relation 
between concentration, mass segregation and core collapse. 
Heating supresses this behavior for Wo = 5. 



2.4 - 



2 - 




2 4 

R 

Figure 5. The radial dependence of the mass spectral index for 
a cluster dominated by bulge shocking: Wo = 5, Mc = IO^Mq, 
M{xp) = 0.8. ai is the initial spectral index, i?^ is the half-mass 
radius. Tidally disrupting clusters may show evidence of mass seg- 
regation. In this case bulge shocking suppresses core contraction, 
leading to expansion and disruption. The profile is approximately 
Wo = 4 and the remaining mass is Mc = 2.3 X 10^ Mq. 



Evaporation times decrease with increasing niu (Figure 
. The decrease in tcv with increasing mass range results 
from enhanced relaxation caused by a more extreme mass 
segregation instability. A 25% range in the duration of strong 
stellar evolution for a = 2.35 gives a range in mass limits 
of 1.9 < ruu < 2.2 and results in very small differences in 
evaporation time. 
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Table 5. Times of core collapse and evaporation 



Wo = 5 




a 






Mc (Mq) 


1.35 


2.35 


3.35 






5.0 X 10^ 


2.7 


3.2 


5.3 


heated 




4.4 


6.3 


14.1 


unheated 


1.0 X 106 


4.3 


5.2 


8.2 


heated 




8.6 


12.3 


25.9 


unheated 


tec 


5.0 X 10^ 


1.6 


1.6 


2.1 


heated 




2.2 


1.9 


2.3 


unheated 


1.0 X 10^ 


2.7 


2.7 


3.6 


heated 




4.0 


3.1 


4.4 


unheated 



Wo = 7 




a 






Mc (Mq) 


1.35 


2.35 


3.35 






5.0 X 10^ 


2.8 


4.0 


7.4 


heated 




4.3 


6.1 


14.0 


unheated 


1.0 X 10'^ 


4.9 


6.4 


11.8 


heated 




8.1 


12.1 


26.8 


unheated 


tec 


5.0 X 10^ 


0.56 


0.37 


0.41 


heated 




0.59 


0.37 


0.38 


unheated 


1.0 X 10'^ 


1.06 


0.71 


0.76 


heated 




1.06 


0.71 


0.76 


unheated 


nil = O.IMq, 


rriu = 


2.OM0 







Table 6. Mass at core collapse 



Wo = 5 




a 






Mc (Mq) 


1.35 


2.35 


3.35 




5.0 X 10^ 


0.39 


0.54 


0.59 


heated 




0.65 


0.88 


0.93 


unheated 


1.0 X 10'^ 


0.30 


0.43 


0.50 


heated 




0.65 


0.88 


0.93 


unheated 



Wo = 7 




a 






Mc (Mq) 


1.35 


2.35 


3.35 




5.0 X 10^ 


0.82 


0.90 


0.92 


heated 




0.93 


0.98 


0.99 


unheated 


1.0 X 10"^ 


0.78 


0.87 


0.90 


heated 




0.93 


0.98 


0.99 


unheated 



mi = O.IMq, niu = 2.0Mq 



5.3 Influence of Anisotropy 

Another internal property that determines the influence of 
external heating is the anisotropy of the stellar orbit dis- 
tribution. Figure ^ shows the variation of heating rate with 
anisotropy radius within a cluster. Heating increases with 
anisotropy due to efficient impulsive heating of radial or- 
bits at apocenter. However, heating rates for ra ~ 2.5 un- 
bind orbits with E > —0.25 in one orbital period tcr and 



quickly alter the DF. The relaxation time is roughly 100 
crossing times, so diffusion cannot maintain the assumed 
level of radial anisotropy in the cluster halo. We estimate 
that anisotropy radii of ra > 5.0 are sustainable through re- 
laxation. The interplay between heating and anisotropy seen 
here provides strong incentive to study the evolution of fully 
anisotropic DFs. 
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Figure 6. Central density evolution for Wo = 5, 10^ clus- 
ters on circular orbits with a = 2.35, mj = 0.1 and niu as in- 
dicated. Evaporation occurs at the termination of each central 
density curve. Evaporation times vary by no more than 10% for 
the expected range 1.9 < m„ < 2.2 in initial upper mass limit 
for a = 2.35. Evaporation times decrease with increases mass ruu 
due to the enhanced mass segregation instability. 



Figure 8. Evaporation times vs. apogalacticon for Wo = 5 (solid) 
and Wo = 7 (dashed). Low mass clusters evaporate within 10 
Gyr in inner Galaxy to apogalactic radii as shown. Strong heating 
drives low eccentricity clusters to evaporation while high densities 
of tidal limitation drive high eccentricity clusters to evaporation. 
Evaporation of k = 0.3, IO^Mq occurs out to average radii of 15 
kpc. 
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Figure 7. Heating rates are shown for clusters on circular orbits 
with indicated anisotropy radius. For ra = 2.5, the anisotropy 
parameter (3 = 1 — Vg/v'^ = 0.17 at the half-mass radius. Energy 
input increases due to efficient impulsive heating of radial stellar 
orbits at apocluster. 



5.4 Evolution in the Milky Way 

5.4-1 Scaled evaporation times 

The dimensionless evaporation times for tidally-limited clus- 



ters discussed in 



5.1 



may be scaled to physical units using 

(16) 



the following relation 

^phys ~ 1-1 X 10 t X tevj 

where i is the orbital period at the tidal radius 



and 



rt 



1/3 



(17) 



(18) 



The quantity np(K,_Ra) is the effective perigalactic angular 
frequency of an orbit of given k and apocentric radius Ra 
due to tidal strain and centrifugal force (defined in Appendix 

§• 

As an example, the dimensionless evaporation times 
given in Table ^ are scaled to a range of apogalactica in 
Figure ^. Clusters evaporate over a wide range of Galacto- 
centric radii depending on In 10 Gyr, clusters on circular 
orbits evaporate within 3 kpc, while those on k = 0.3 orbits 
evaporate out to average radii of 15 kpc. 



5.4.2 Survival and disruption 

Here we present a simple evolutionary scenario in which clus- 
ters form at apocenter with a range of mean density param- 
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eterized by Pcrit.FR{R), the Fall-Rees (1985) critical cloud 
density at Galactocentric radius R. This parameterization 
is chosen to allow normalization with respect to a particular 
model. Other models (e.g. Harris & Pudritz 1994; Murray 
& Lin 1992) can be similarly evaluated given expressions 
for initial protocloud densities as a function of Galactocen- 
tic radius. A range of density is used to define a range of 
M{xp) for clusters at each radius, thereby illustrating char- 
acteristics which are independent of any particular model. 
We only consider relaxation, external heating and binary 
heating although gas removal and stellar evolution will play 
an important role following formation. These effects should 
weaken the potential and increase disruptive tendencies de- 
scribed here. 

In the first case, clusters form on eccentric k = 0.3 orbits 
(e.g. Eggen, Lynden-Bell and Sandage 1962). Figure ^ shows 
the resulting pattern of survival, disruption and evaporation 
for 10^ and IO^Mq clusters after 10 Gyr. Clusters initially 
with 10^ Mq do not surive within Rav = 15, reflecting the 
evaporation times shown above. Lower density clusters suffer 
disruption to even larger distances. High mass clusters with 
M{xp) < 0.8 can suffer disruption but none can evaporate. 

Cluster formation on less eccentric n — 0.7 orbits shows 
the same qualitative pattern of survival, evaporation and 
disruption as above (Figure ^o|). The consequences are less 
severe because the density contrast between formation at 
apocenter and tidal limitation at pericenter is not as great. 
In this case, low mass cluster survival is limited to regions 
beyond 5 kpc for clusters which are nearly tidally-limited. 



6 IMPLICATIONS FOR MILKY WAY 
CLUSTERS 

The calculations presented above bear on our understand- 
ing of the observed mass and space distributions of Galactic 
globular clusters. We summarize some relevant properties 
for reference. In the Djorgovski (1993) compilation, 65% of 
130 clusters with distance estimates lie within the solar cir- 
cle. The overall peak of the luminosity function of Galactic 
globulars is M„ = —7.36 (Harris 1991) corresponding to a 
mass of 1.5 x IO^Mq (where M/L„ = 2). The luminosity 
function varies little in this inner region. 

Our results imply that the observed characteristics 
of this inner population have evolved with time. Because 
10^ Mq clusters evaporate or lose large amounts of mass in 
a Hubble time in the inner Galaxy, clusters at the peak of the 
luminosity function have evolved from higher mass. For ex- 
ample, at 6 kpc clusters on circular orbits with Mv — —7.36 
will evolve to AIv = —6.8 in 10 gyr, losing roughly 40% of 
their initial mass. Inside the solar circle, clusters near the 
present peak had at least 30% more mass, depending on the 
orbit. 

Many clusters will also have vanished. We predict that 
evaporation and disruption of 10^ M© clusters occur within 
3 kpc for K — 1.0 and within apocentric radii Ra ~ 20 kpc 
for K — 0.3. For intermediate k, the destruction region is 
bracketed by these limiting cases. These results buttress ar- 
guments based on two-body relaxation times that the shape 
of the luminosity function stems from evaporation and dis- 
ruption of a larger initial population of low mass clusters 
(e.g. Larson 1996, Okazaki & Tosa 1995). 



The dependence of survival on orbit implies that the 
kinematic distribution of clusters has evolved as well. Clus- 
ters on high eccentricity orbits in the inner Galaxy are least 
likely to survive due to both evaporation and bulge shock- 
ing. This suggests a decrease in the net velocity dispersion 
for the rotating system of metal-rich and inner halo metal- 
poor clusters (Zinn 1993). This tendency may also partially 
account for observed differences between the kinematics of 
halo field stars and metal-poor globular clusters (e.g. Aguilar 
et al. 1988). 

Finally, survival also depends strongly on initial cluster 
densities. Destruction is more pronounced for clusters with 
low initial density and low initial concentration due to bulge 
shocking. Bulge shocking can disrupt massive 10^ Mq clus- 
ters on eccentric orbits out to 40 kpc provided M{xp) < 0.8. 
However, a proper assessment of the initial distribution of 
cluster densities requires cosmogonical considerations. 

We conclude that the segment of the cluster population 
which is native to the Milky Way or which was accreted 
at an early time represents a dynamically selected sample, 
with current masses, orbits and densities all favored for sur- 
vival over a Hubble time of evolution. Tidal interaction with 
the Galactic disk will amplify these effects. Details will be 
described in a subsequent paper. 



7 SUMMARY 

The key conclusions are as follows: 

(i) Time-dependent heating on low-eccentricity orbits ac- 
celerates evolution and sharply reduces evaporation times. 

(ii) Tidally limited clusters on high-eccentricity orbits 
have high internal density, leading to short evaporation 
times even though heating rates are negligible. 

(iii) Bulge shocking on high-eccentricity orbits can 
rapidly disrupt clusters over a wide range in mass and 
apogalactic radius when their densities are roughly a factor 
of 10 below the mean density required for tidal limitation. 

(iv) Evaporation and disruption have shaped the mass, 
orbit and density distribution of clusters. In particular, clus- 
ters at the peak of the luminosity function had at least 
~ 30% more mass depending on orbit. Evaporation on high- 
eccentricity orbits has decreased the velocity dispersion in 
the cluster kinematic distribution. 

Secondary results are as follows: 

(i) Evaporation times do not strongly depend on concen- 
tration in most cases. However, heating can lead to rapid 
disruption in massive clusters with low concentration be- 
cause of the low binding energy and long relaxation time. 

(ii) Clusters disrupting due to heating may still show 
signs of mass segregation due to continued relaxation. 

(iii) Heating accelerates evolution over a range of mass 
spectral index and reduces the dependence of evaporation 
time on different initial mass spectra. 

(iv) The development of anisotropy through relaxation 
in the core will increase evolutionary rates found in the 
isotropic distributions investigated here. 
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Figure 9. Contours of log mass (solid) show survival and disruption of proto-clusters with an initial Wo = 5 profile after 10 Gyr due 
to tidal effects onK = 0.3 orbits. Contours arc in the range 3.75 < logM < 4.5. Rav indicates average orbital radius and Pc-u FB. is the 
Fall-Rees (1985) critical cloud density at radius R. Dotted contours show lines of equal M{xp). Left: 10^ Mq clusters with p Pcrit fr{R) 
evaporate due to high initial densities, lower density clusters disrupt from bulge shocking and clusters with M(xp) ^ 0.95 survive longest 
due to balance between heating and relaxation. Right: a density of 0.1pcrit,FR leads to bulge shocking and disruption in 10^ Mq out to 
Rav ~ lOfcpc. 




Figure 10. As in Fig. 8 but on k = 0.7 orbits. Left: high density IQ^Mq clusters evaporate, low density clusters disrupt due to extreme 
tidal heating and clusters survive at larger radii. The convergence of contours into the upper left corner is a numerical artifact caused 
by mean densities beyond the range of our calculations. However, these clusters also evaporate because of the high densities. Right: low 
density clusters disrupt at densities roughly 30% of the mean density for tidal limitation out to 15 kpc due to strong tidal heating 



The Effect of the Galactic Spheroid on Globular Cluster Evolution 13 



ACKNOWLEDGEMENTS 

This work was supported in part by NASA award NAGW- 
2224. 



REFERENCES 

Aguilar, L., Hut, P. & Ostriker, J. P. 1988, ApJ, 335, 720 
Binney, J. J. & Tremaine, S. D. 1987, Galactic Dynamics, (Prince- 
ton: Princeton U. Press) 
Bracewell, R. N. 1986, Tlic Fourier Transform and its Applica- 
tions, (New York: McGraw Hill) 
Chandrasekliar, S. 1942, Principles of Stellar Dynamics, (Chicago: 

Univ. Chicago Press) 
Chernoff, D. F. & Djorgovski, S. C. 1989, ApJ, 339, 904 
Chernoff, D. F., Kochanek, C. S. & Shapiro, S. L. 1986, ApJ, 309, 
183 

Chernoff, D. F. & Shapiro, S. L. 1987, ApJ, 322, 113 
Chernoff, D. F. & Weinberg, M. D. 1990, ApJ, 351, 121 
Cohn, H. N. 1979, ApJ, 234, 1036 
Cohn, H. N. 1980, ApJ, 242, 765 

Cohn, H. N. 1985, in Dynamics of Star Clusters, eds. Goodman, 

J., and Hut, P., (Boston: D. Reidel Publishing) 
Djorgovski, S., Piotto, G. & Capaciolli, M. 1993, AJ, 105, 2148 
Djorgovski, S., & Meylan, G. 1994, AJ, 108, 1292 
Drukier, G. A., Fahlman, G.G. & Richer, H. B. 1992, ApJ, 386, 
106 

Eggen, O. J., Lynden-Bell, D. & Sandage, A. 1962, ApJ, 136, 748 
Elson, R. A. W., Fall, S. M. & Freeman, K. C. 1987, ApJ, 323, 

54 

Elson, R., Hut, P. & Inagaki, S. 1987, ARAA, 125, 565 
Fall, S. & Rees, M. 1985, ApJ, 298, 18 

Goldstein, H. 1980, Classical Mechanics, (Reading, Mass: 

Addison- Wesley), p.383 
Harris, W. E. 1991, ARAA, 29, 543 
Harris, W. E. & Pudritz, R. E. 1994, ApJ, 429, 177 
Hcggie, D. C, & Ramamani, N., 1995, MNRAS, 272, 317 
Hcnon, M., 1961, Ann. d'Astrophys., 24, 369 
Hernquist, L. & Ostriker, J. P. 1992, ApJ, 386, 375 
Hcsser, J. E., 1993, in The Globular-Cluster Galaxy Connection, 

eds. Smith, G. H. & Brodie, J. P., (Astr. Soc. Pacific: San 

Francisco) 

Inagaki, S., in Dynamics of Star Clusters, eds. Goodman, J., and 

Hut, P., (Boston: D. Reidel PubUshing) 
Inagaki, S & Saslaw, W. 1985, ApJ, 292, 339 
King, I. R., 1962, AJ, 67, 471 

Krall, N. A., & Trivelpiece, A. W. 1973, Principles of Plasma 

Physics, (New York: McGraw-Hill) 
Landau, L. D., & Lifschitz, E. M. 1965, Quantum mechanics, 

non-relativistic theory, (London: Pergamon Press) 
Larson, R. B. 1990, PASP, 102, 709 
Larson, R. B. 1994, preprint 
Lee, H. M. & Goodman, J. 1995, ApJ, 443, 109 
Lee, H. M. & Ostriker, J. P. 1987, ApJ, 322, 123 
Lee, H. M., Fahlman, G. G., & Richer, H. B. 1991, ApJ, 366, 455 
Lichtcnberg, A., & Lieberman, M. 1983, Regular and Stochastic 

Motion, (New York: Springer- Verlag) 
Lightman, A. & Shapiro, S. 1978, Revs. Mod. Phys., 50, 437 
Lynden-Bell, D. & Kalnajs, A. 1972, MNRAS, 157, 1 
Lynden-Bell, D. & Wood, R. 1968, MNRAS, 138, 495 
Murray, S. D. & Lin, D. N. C. 1992, ApJ, 400, 265 
Oh, K. S. & Lin, D. N. C. 1992, ApJ, 386, 519 
Okazaki, T. & Tosa, M. 1995, MNRAS, 274, 48 
Ostriker, J. P.1985, in Dynamics of Star Clusters, eds. Goodman, 

J., and Hut, P., (Boston: D. Reidel Publishing) 
Ostriker, J. P., Spitzer, L., Jr. & Chevalier, R. 1972, ApJL, 176, 

L51 



Richer, H. B., & Fahlman, G. G. 1989, ApJ, 339, 178 
Rosenbluth, M. N., MacDonald, W. M. & Judd, D. L. 1957, Phys- 
Rev, 107, 1 

Spitzer, L., Jr. 1987, Dynamical Evolution of Globular Clusters, 

(Princeton: Princeton Univ. Press) 
Spitzer, L., Jr. & Hart, M. 1971, ApJ, 164, 399 
Spitzer, L., Jr. & Chevalier, R. A. 1973, ApJ, 183, 565 
Stix, T. H. 1992, Waves in Plasmas, (New York: AIP Press) 
Szebehely, V. 1967, Theory of Orbits, (New York: Academic 

Press) 

Titchmarsh, E. C. 1986, Introduction to the Theory of Fourier 

Integrals, (New York: Chelsea Publishing Company) 
Tremaine, S. D. & Weinberg, M. W. 1984, MNRAS, 209, 729 
Weinberg, M. D. 1994, AJ, 108, 1414 

Zinn, R., 1993, in The Globular-Cluster Galaxy Connection, eds. 
Smith, G. H. & Brodie, J. P., (Astr. Soc. Pacific: San Fran- 
cisco) 



APPENDIX A: DERIVATION OF TIDAL 
POTENTIAL 

In the inertial Galactocentric frame, the coordinate compo- 
nents of a cluster star are 



R — ^ ~t~ Rcomi 

and its velocity components are 

V = V + Vcom, 



(Al) 



(A2) 



where r and v are the coordinates and velocities of the mem- 
ber star measured relative to the cluster center of mass and 
Room and Vcom are the ccnter-of-mass position and veloc- 
ity of the cluster. The Hamiltonian for an individual star in 
these coordinates is therefore 



Ho = ^\Vf + M\R - Rcom.\) + '^r{\R\) 



(A3) 



Wc introduce coordinates centered on the cluster with 
axes fixed in space through a canonical transformation using 
a generating function of the second kind (Goldstein 1985). 
This function can be written 



F2{R, V, t) = {v + Vcom) ■ {R - Rcom) + /(*) 



(A4) 



where fit) is an arbitrary function of time. The transforma- 
tion obeys the conditions Vi = dF2/dRi and = dF2/vi, 
thus satisfying Hamilton's principle. The new Hamiltonian 
Ho = Hq + dF2 /dt so that (assuming the summation con- 
vention throughout) 



^^0 = ^1^1^ + *c(|rl) + ^G{\f+ Rcom\) 



dRi 



-\V P + ^ 

^\Vcom\ + g^. 



(A5) 



Expanding the Galactic tidal potential about the center of 

mass, we obtain 



Ho = hv\' + M\n)+[ 



1 d^^G 



2 dRiR, 



Tir-j + ...] 



Rao 



+ [-\\Vcom\''+^G{\Rcom\)] + ^. (A6) 
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The term in the second pair of brackets is an arbitrary func- 
tion of time which arises as an ambiguity in canonical trans- 
formations (Goldstein 1985). We note that it equals —Lcom, 
the negative Lagrangian of the center-of-mass motion and 
that it can be eliminated by an appropriate choice of /. In 
this case, setting / = / Ldt (the action associated with the 
center-of-mass motion) gives the desired form of the Hamil- 
tonian for a star in the cluster frame: 



Ho = llvY 



+ <i>e(|r1) + [i 



2 dR^Ri 



(A7) 



In the expansion of the tidal potential, we can ignore 
all but the lowest-order term because successive terms are 
proportional to (|r|/|-Rcom|)"~'^ {n = 3,...) relative to the 
second order term and, therefore, fall off quickly due to 
the small size of a cluster compared to the size scale of the 
Galaxy. We are thus left with the quadratic approximation 
to the Galactic tidal field. The expression for the perturbing 
potential, equation (8) , is obtained by evaluating this term of 
the expansion for the specific case of the logarithmic sphere. 



APPENDIX B: PERICENTRIC INNER 
LAGRANGE POINTS 

To obtain an expression for the pericentric inner Lagrange 
point, it is convenient to first transform to rotating coordi- 
nates having one axis aligned with the galactocentric radius 
of the cluster on its orbit. We omit the details of the transfor- 
mation here and simply give the expression for the effective 
potential at pericenter in the rotating frame: 



$,(|r1) + if7g(i?,)[|rf -2:r^T - 



(Bl) 



The first term is the cluster potential. The second term is the 
quadratic tidal potential for the logarithmic sphere trans- 
formed to a rotating coordinate system. The last term is 
the centrifugal potential arising from the angular frequency 
of rotation at pericenter. The quantity f2p is the pericentric 
angular frequency of the cluster orbit while f2o(^?p) is the 
angular frequency of a circular orbit at the pericentric radius 
which defines the tidal strain. 

The pericentric inner Lagrange point Xp occurs at 
the instantaneous inflection point in the effective potential 
which lies along the galactocentric radius of the cluster. We 
derive an expression for by considering the balance of 
forces which is implied by the effective potential. Taking the 
gradient, and considering the instantaneous point of equi- 
librium along the Galactocentric radius gives 



(B2) 



where the left-hand side gives the cluster force while the 
right-hand side gives the centrifugal force and tidal strain, 
respectively. 

To derive an expression for the inner Lagrange point in 
terms of Ra/Rp, we first rewrite the angular frequencies Q,p 
and Q.o{Rp) in terms of the angular frequency of a circular 
orbit at apocenter Q.o{Ra). Using conservation of angular 
momentum between apocenter and pericenter gives 



J I Rp — 



V^oiRa] 



Ra 
Rv 



(B3) 



where r} denotes the ratio of the angular frequency of the 
orbit to the angular frequency of a circular orbit with the 
same apocenter. The flat rotation curve defines 



no{Rp) = ^ = no{Ra){^) 



and substituting into equation 
obtain 



(B4) 

and solving for Xp, we 



GM{Xp) \l ( Ra 



2Q.l{Ra) \2\R, 



(B5) 



Now we define the effective pericentric angular frequency 



n'^^nl{Ra){\(^^ 



2 I Ra , , , 



(B6) 



where k = rie^^ ^ ''^^ for the logarithmic sphere. 



APPENDIX C: DERIVATION OF FLUX 
EQUATION 

The flux equation, equation can be derived from the (/2) 
formalism given in Weinberg (1994). The function (/2) de- 
fines the externally induced change in the distribution func- 
tion. The general form for {/2) is 



(Cl) 



where we have performed the phase-averaging to derive the 
perturbation as a function of the actions. The quantity Wi(I) 
is a scalar function of the actions 



M^.(I) = i(l-§)lU.na(l.fl)^ 

where a(I ■ 12) is a Fourier coefficient given by 



?(t') 



(C2) 



(C3) 



a(I ■n)= dt'e" 

J to 

and other quantities are as defined in §^ 

To derive the fiux equation we note that {/2) can be 
written as a divergence: 



;/2) = E^-W' 



(C4) 



where Wi = Wi x 1. This equation makes number conserva- 
tion manifest in action space. 

To implement this term in a 1-dimensional Fokker- 
Planck scheme, we must change variables from actions to 
{E, K, cos P) and average over k and cos /3 to obtain the one- 
dimensional fiux in energy space. The transformation can 
be performed easily using the covariant form of the equa- 
tion (e.g. Rosenbluth et al 1957). 

The divergence written in covariant form is 



_i d_ 



(C5) 



where is a contravariant vector and g is the determinant 
of the metric tensor, equal to the square of the Jacobian. 
Using this equation, we transform to the new [E, k, cos 13) 
coordinates, which we denote using primes. Transforming 
the contravariant vector Wi gives 
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wi = Wi(i- n)E + Wi 



I2 



Jmax(E) Q,\Q,2i\ 
h 



+Wi 



COS /3, 



(C6) 



where the quantity Wi is the function defined above but now 
written in terms of the new variables. Wj is the function 
equivalent to Wi in the new coordinates. 
Noting that the Jacobian is 



we may write {/2) in the new coordinates: 



max 



Averaging over (k, cos/3) 

_ / dndcos /3Kj^ax/^i{/2) 

J dud cos PK,J^ax/^l 

gives the total change as a function of energy: 

/ dKdcOSf3Kjla:c/^ld{j2Wlil- n))/dE 



J dud cos /3KJmax/f2l 



(C7) 



(C8) 



(C9) 



(CIO) 



where the fluxes in the k and cos /3 directions vanish due to 
the averaging. Fully expressed, the equation reads 

- { ^ E ^ [/ ■§)('• «)iuii>(i ■ 

i^J?nax/^i]^^J dnnj^^^/ni^ (Cll) 
The phase space volume 



16tt^P{E) = (27r) 



r2 

3 / ^^max I 



(C12) 



which we substitute to find the total phase space flux 
((/2)) = im^P(E)r'^{{{E)}}. (C13) 

In the asymptotic limit, this becomes the rate of change of 
the phase space density, equation (ph. 



APPENDIX D: IMPLEMENTATION 

The rate of change in the distribution function due to ex- 
ternal heating is given by equation (^. We write this in fi- 
nite difference form for consistency with the Fokker-Planck 
scheme and solve after each diffusion step. The numeri- 
cal implementation uses a flux- conserving finite-difference 
scheme with explicit time advance: 



J.i 



1 



At 



•^7+1/2 — ^■ 



AE 



where the flux is denoted 
T -R '^f 



(Dl) 



(D2) 



0.015 




Acj = 0.0025 



-0.4 -0.3 -0.2 -0.1 
E 



0.015 ■ 



0.01 



0.005 - 



Acj = 0.005 



-0.4 -0.3 -0,2 -0.1 
E 




Figure Dl. Increasing grid spacing corresponds to broadening 
resonances for the resonant heating calculation. Broader reso- 
nances spread the input power over a range of energies. As a 
result, the DF evolves more slowly and does not develop strong 
resonant peaks and troughs. Here we evolve the DF in a fixed 
potential 



and we have rewritten the previously defined heating rate 



{{£)) = R{E) 



dl_ 
dE' 



(D3) 



The function R{E) represents the sum over all resonances 
which couple to orbits of energy E. 

In equation (^, 5- functions denote resonant coupling 
of internal and external orbital frequencies. However, the 
resonant interaction has finite duration due to non-linear 
saturation or detuning which corresponds to a width in fre- 
quency space. For weak perturbations, narrow resonances 
develop since orbital frequencies evolve slowly. For strong 
perturbations large widths occur because frequencies evolve 
rapidly. 

The grid spacing employed in the difference scheme de- 
fines the frequency widths. Wider spacing implies broader 
resonances. Broader resonances reduce the heating rate by 
smearing the input power over a wide region in phase space 
as shown in Figure 



Dl 



To estimate the proper grid spacing, 
we use the bandwidth theorem or uncertainty relation (e.g 
Bracewell 1986). The bandwidth theorem identifies the re- 
ciprocal relationship between the frequency width and sam- 
pling time of an oscillator: 

Alu = =^. (D4) 
At ^ ' 

Here the sampling time is the duration of resonance. Since 
non-linearity develops with some typical change in energy 
&E, we use the rate of energy input to estimate the duration 
of resonance: 

= ^. (D5) 
E 
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Comparisons with simulation indicate that 10% change in 
energy typically leads to frequency evolution. We calibrate 
the appropriate frequency widths using fully self-consistent 
N-body calculations described in Appendix ^. For typical 
heating rates in tidally limited clusters, an initial spacing of 
Au! = AI- f2 ~ 0.005 is appropriate. We use this value for all 
tidally-limited calculations. Larger tidal truncations require 
an increased width. 

Two final implementation issues are the boundary con- 
ditions on the external heating equation and the Fokker- 
Planck equation. For the boundary conditions on equation 
(Dl), we set the flux to zero at the center and the gradient 
of the flux to zero at the edge. The latter condition repre- 
sents evaporation. The last grid point of the DF therefore 
stays fixed between each diffusion step. We tested the choice 
in outer boundary condition using a zero flux condition and 
found solutions which differed by a few percent at most. 

We use the standard tidal boundary condition in the 
Fokker-Planck equation. This calls for truncating the distri- 
bution function at the maximum energy allowed by the tidal 
limit. Because strict application of the boundary condition 
calls for truncating the cluster at Xp, we would throw out 
a potentially large fraction of the initial cluster distribution 
in cases where M{xp) < 1. As a less extreme alternative, we 
place the zero-DF boundary at the initial limiting energy of 
the cluster and set the heating rate according to our choice 
of Xp. Then we allow the boundary to evolve according to 
the mass loss. Tests with N-body calculations shown in Ap- 
pendix |e| indicate that the prescription works correctly. 
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Figure El. A Wo = 5 cluster on a circular orbit with M{xp) = 
0.95. The half-mass dynamical time t^y„ = 0.3 and the orbital 
period is 2it. The top left panel shows the total mass as a function 
of time while the remaining panels show the mass profile at the 
indicated times. Agreement is good for lOOt^yn ■ Deviation at later 
times results from inherent non-linearity and relaxation. 



APPENDIX E: COMPARISON WITH 
SIMULATION 

Comparisons of the external heating theory with N-body 
simulations are used to test linearity, the assumption of 
isotropy over long times, the importance of non-spherical 
moments in the potential, and t he b oundary condition on 



the heating equation (equation Dl) and to calibrate fre- 
quency widths of heating rates. We use a self-consistent 
field expansion code (e.g. Hernquist & Ostriker 1992) with 
1.5 X 10"* particles, radial expansion order n = 10 and an- 
gular order 1 = 4. Figures ^ and E2 show comparisons for 



K = 1.0 and K = 0.3 orbits. The two methods agree fairly 
well, especially at early times, before non-linearity and re- 
laxation produce differences at later times. 
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Figure E2. The same cluster as above on a re = 0.3 orbit with M{xp) = 0.95. The orbital period is t = 30. The left panel shows the 
evolution in total mass while the right panel compares the mass profiles at t = 145 or 483tdj,„. Agreement is good over this duration. 



